In silico design and immunoinformatics analysis of a universal multi-epitope vaccine against monkeypox virus

Monkeypox virus (MPXV) outbreaks have been reported in various countries worldwide; however, there is no specific vaccine against MPXV. In this study, therefore, we employed computational approaches to design a multi-epitope vaccine against MPXV. Initially, cytotoxic T lymphocyte (CTL), helper T lymphocyte (HTL), linear B lymphocytes (LBL) epitopes were predicted from the cell surface-binding protein and envelope protein A28 homolog, both of which play essential roles in MPXV pathogenesis. All of the predicted epitopes were evaluated using key parameters. A total of 7 CTL, 4 HTL, and 5 LBL epitopes were chosen and combined with appropriate linkers and adjuvant to construct a multi-epitope vaccine. The CTL and HTL epitopes of the vaccine construct cover 95.57% of the worldwide population. The designed vaccine construct was found to be highly antigenic, non-allergenic, soluble, and to have acceptable physicochemical properties. The 3D structure of the vaccine and its potential interaction with Toll-Like receptor-4 (TLR4) were predicted. Molecular dynamics (MD) simulation confirmed the vaccine’s high stability in complex with TLR4. Finally, codon adaptation and in silico cloning confirmed the high expression rate of the vaccine constructs in strain K12 of Escherichia coli (E. coli). These findings are very encouraging; however, in vitro and animal studies are needed to ensure the potency and safety of this vaccine candidate.

Introduction After more than two years of significant global economic and health impacts from COVID-19, outbreaks of the monkeypox virus (MPXV) have been reported in many countries around the world, raising new concerns about a new global pandemic [1]. World Health Organization (WHO) has reported 87,301 confirmed cases and 130 confirmed deaths of MPX in 111 countries as of May 3, 2023 continents [2]. The virus was named MPXV after it was discovered in monkeys in a Danish laboratory in 1985 [3]. The first human case was identified in a ninemonth-old child from the Democratic Republic of the Congo (DRC) in 1970 [4]. MPXV has since become endemic in the DRC and has spread to other African countries. The first MPXV outbreak outside of Africa has been reported in the United States, following the shipment of infected African mammals from Ghana [5].
The MPXV is divided into two clades: the Congo Basin (central African) clade and the West African clade. The Congo Basin clade is more virulent, with a fatality rate of 106% compared to the West African clade, which has a fatality rate of 3.6% [6]. MPXV is a zoonotic virus that can be transmitted from animal to human and from human to human [7]. Animal-to-human transmission can occur through biting or scratching, consumption of game meat, and direct or indirect contact with body fluids or feces [8]. Transmission from human to human may occur via large respiratory droplets. Because respiratory droplets can only travel a few feet, prolonged face-to-face contact is required to spread the virus. MPXV can be transmitted through contact with an infected person's body fluids, feces, and urine, as well as infected clothing or bed linen [1]. The MPXV virus may be released in wastewater during rinsing, showering, or using the toilet [9]. The current outbreak of this virus among many gay men or men who have sex with other men raises the possibility of sexual transmission of the virus [10].
The symptoms of MPXV infection are similar to but milder than those of smallpox [11] and include fever, chills, headache, muscle aches, and fatigue. The main difference between this disease's symptoms and those of smallpox is that MPX causes lymph node swelling (lymphadenopathy) [1]. The incubation period of MPX is usually 7-14 days, but it can take up to 21 days [1]. The infected person develops skin rashes that resemble blisters 1 to 3 days (sometimes more) after the onset of fever. The majority of these rashes begin on the face and spread to other areas of the body [12]. MPXV infection is diagnosed based on the patient's history, clinical symptoms, and laboratory tests such as PCR, ELISA, western blot, and immunohistochemistry. The World Health Organization (WHO) recommends using qRT-PCR to detect viral DNA during acute infection of the MPXV [1].
MPXV belongs to the Poxviridae family's Orthopoxvirus genus [13]. The morphology of the MPXV reveals that virions are ovoid or brick-shaped particles enclosed by a geometrically corrugated lipoprotein outer membrane, as with other orthopoxviruses. The size range of the MPXV is 200 by 250 nm [14]. The MPXV genome is a linear double-stranded DNA with a length of about 197 kbp, consisting of hairpin loops, some open reading frames, and tandem repeats, and both ends contain inverted terminal repeats (ITRs) with a length of 10 kbp [15].
Vaccines have historically been the most effective tool for preventing and even eradicating infectious diseases. Vaccines can also cause herd effects, which result in protection even among those who have not been vaccinated, which can be especially useful for low-income individuals who are unable to access health care [16]. So far, no vaccine has been developed specifically to protect against MPX infection. The FDA approved ACAM2000 for use against smallpox in 2007. Since this vaccine is made from a live, replication-competent Vaccinia virus, there is a risk of serious side effects in vaccinated individuals [17]. On the other hand, Jynneos, a non-replicating modified Vaccinia Ankara virus vaccine, was approved by the FDA in 2019 for the prevention of both MPX and smallpox. Jynneos does not cause the production of live virus in vaccinated individuals, making it safer than ACAM2000, particularly for use in immunocompromised individuals; however, there is limited information on Jynneos's effectiveness in preventing MPXV infection in humans [18].
The emergence of new infectious diseases necessitates the development of novel vaccine design methods. Conventional vaccine development methods are time-consuming and costly, requiring the cultivation of pathogenic microorganisms and the identification of their immunogenic components [19]. Bioinformatics methods and techniques have expanded exponentially over the years, recently influencing immunological studies. Reverse vaccinology is a computational approach that uses bioinformatics methods to identify new antigens from genome sequence information without the need to isolate and culture the pathogen, significantly accelerating vaccine development [20]. The design of multi-epitope vaccines using reverse vaccinology is a new and rapidly growing field. These vaccines have gotten a lot of attention because of their high specificity, safety, low-cost production, and simultaneous induction of cellular and humoral immune responses [21,22]; the major drawback of these vaccines is their low immunogenicity [23][24][25], which is why it is suggested that protein-based adjuvants be included in the structure of these vaccines to solve this problem [21].
This study aims to design a multi-epitope vaccine against the MPXV using a reverse vaccinology approach. The cell surface-binding protein and envelope protein A28 homolog were used as target proteins for epitope prediction in this study. The cell surface-binding protein and envelope protein A28 homolog are attractive targets for vaccine development for the following reasons: the cell surface-binding protein is responsible for important biological functions of the MPXV, including virion attachment to the host cell and entry of the virus into the host cell [26]; the envelope protein A28 homolog is also required for virus entry into the host cell and cell-cell fusion [27]. To improve the vaccine candidate's immunogenicity, the cholera toxin B subunit (CTxB) sequence, which has been proven to act as a potential viral adjuvant, was incorporated into the vaccine construct [28][29][30]. We hope that the findings of this study will help in the fight against MPXV by leading to the development of a vaccine.

Materials and methods
The sequential steps used in this study to design and evaluate a multi-epitope vaccine against the MPXV are depicted in Fig 1.

Retrieval of protein sequences and identification of conserved regions
For each cell surface-binding protein and envelope protein A28 homolog, a hundred sequences in FASTA format were obtained from the NCBI database (https://www.ncbi.nlm. nih.gov/). Following the removal of irrelevant sequences, partial sequences, and sequences with ambiguous characters, multiple sequence alignment (MSA) was done using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) with default parameters [31].

Prediction and screening of T-cell and linear B-cell epitopes
Cytotoxic T lymphocyte (CTL) epitopes play an important role in the elimination of virusinfected cells, the promotion of cellular immunity, and the reduction of virus levels in the body [32]. The ProPred-I server (https://webs.iiitd.edu.in/raghava/propred1/index.html) was used to predict CTL epitopes from the cell surface-binding protein and envelope protein A28 homolog of MPXV. ProPred-I is an online web tool that predicts MHC binding sites in an antigenic sequence for 47 MHC class-I alleles using a matrix-based approach [33]. The adaptive immune response relies heavily on helper T-lymphocytes. They play roles in the activation of B-cells and the killing of infected target cells, respectively [34]. The ProPred server (https:// webs.iiitd.edu.in/raghava/propred/index.html) was used to predict the helper T lymphocyte (HTL) epitopes of the target proteins. ProPred predicts HTL epitopes in an antigen sequence using quantitative matrices for 51 HLA-DR alleles, which cover more than 90% of MHC Class II molecules expressed on antigen-presenting cells [35]. During a viral infection, B-cells use The process comprises numerous key phases, including target protein sequence retrieval, epitope prediction, vaccine construction, physicochemical characterization, 3D model prediction, molecular docking to immune receptor, MD simulation, and in silico cloning.
The vaccine components, while antigenic and non-toxic, should not cause allergic reactions. As a result, the predicted epitopes were tested for antigenicity, toxicity, and allergenicity. VaxiJen 2.0 server (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) was used to calculate the antigenicity of the predicted epitopes. This is the first server for independent prediction of the level of protective antigens, allowing antigen classification based solely on the physicochemical properties of proteins without the use of sequence alignment [40][41][42]. In this study, we chose the virus as the target organism and set the antigenicity threshold at 0.4. The ToxinPred server (https://webs.iiitd.edu.in/raghava/toxinpred/design.php) was used to predict the potential toxicity of epitopes. A key feature of the server is that it calculates various physicochemical properties in addition to toxicity [43,44]. The AllerTOP v. 2.0 server (https:// www.ddg-pharmfac.net/AllerTOP/) was used to predict the allergenicity of the epitopes. To predict allergenicity, this server employs the auto cross covariance (ACC) transformation of amino acid sequences into uniform vectors of equal length [45]. The presence of predicted epitopes in the target proteins' conserved regions was then investigated. HTL epitopes that induce cytokines are important for vaccine development [46]. In addition to the previously mentioned parameters, HTL epitopes were tested for the production of cytokines such as interferongamma (IFN-γ) and interleukin-4 (IL-4). The two cytokines, IL-4 and IFN, play important roles in the generation and regulation of immune responses. IFN-γ promotes T helper type 1 responses, while IL-4 promotes T helper type 2 responses [47]. The IFNepitope server (https:// webs.iiitd.edu.in/raghava/ifnepitope/design.php) was used to predict the ability of HTL epitopes to induce IFN-γ production [48]. In this study, we chose an SVM-based approach and an IFN-gamma versus other cytokine model for prediction. Furthermore, the IL4pred server (https://webs.iiitd.edu.in/raghava/il4pred/design.php) at 0.2 threshold was used to check HTL epitopes for IL-4 production [49].

Multi-epitope vaccine construction
The final selected CTL, HTL, and LBL epitopes from the cell surface-binding protein and envelope protein A28 homolog were linked to each other using appropriate linkers to design a multi-epitope vaccine construct. The CTL, HTL, and LBL epitopes were connected using AAY, GGGGS, and KK linkers, respectively. Furthermore, to enhance the immunogenicity of the designed vaccine, the CTxB (accession number: P01556) was added as an adjuvant to the N-terminal of the vaccine construct via an EAAAK linker.

Population coverage analysis
The frequency of various HLA alleles varies by geography and ethnic background [50,51]. The IEDB population coverage analysis tool (http://tools.iedb.org/population/) was used to assess the population coverage of the vaccine candidate [52]. We used the selected CTL and HTL epitopes, as well as their MHC alleles, separately and in combination for this purpose. We additionally emphasized the total coverage of selected alleles across multiple continents.

Prediction of antigenicity, allergenicity, solubility, and physicochemical properties of the vaccine construct
The antigenicity of the vaccine construct was predicted by the VaxiJen 2.0 server (threshold value of 0.4) and the ANTIGENpro server (http://scratch.proteomics.ics.uci.edu/). ANTIGENpro predictions are sequence-based, alignment-free, and pathogen-independent. To make predictions, this server employs a two-stage architecture based on multiple representations of the primary sequence and five machine learning algorithms [53]. The AllerTOP v. 2.0 server was used to predict the allergenicity of the proposed vaccine. To test the solubility of the multi-epitope vaccine, the Protein-Sol server (https://protein-sol.manchester.ac.uk/) was used. The population average for the experimental dataset in this server is 0.45, so if the solubility value of a protein sequence is predicted to be greater than 0.45, it indicates that it has a higher solubility than the average soluble protein of Escherichia coli (E. coli) from the experimental solubility dataset [54]. Furthermore, the Expasy server's ProtParam tool (https://web. expasy.org/protparam/) was used to predict various physicochemical properties of the final vaccine construct, including the number of amino acids, molecular weight, theoretical pI, halflife, instability index, aliphatic index, and grand average of hydropathicity (GRAVY) [55].

3D structure modeling, refinement, and validation of the vaccine candidate
The I-TASSER server (https://zhanggroup.org/I-TASSER/) was used to predict the threedimensional structure of the proposed vaccine. This server is a web-based resource for predicting protein structure and annotating structure-based functions. I-TASSER first recognizes structural templates from the PDB using multiple threading alignment approaches. Iterative fragment assembly simulations are then used to build full-length structural models. Finally, functional insights are obtained by matching predicted structure models with known proteins in function databases [56][57][58]. The predicted 3D model was refined to improve structural quality using the GalaxyRefine server (https://galaxy.seoklab.org/cgi-bin/submit.cgi?type= REFINEO). CASP10 has successfully tested the GalaxyRefine server's refining approach. The method first reconstructs side chains, then performs side chain repacking and overall structure relaxation using molecular dynamics simulation [59]. Validation is an important step in the modeling of proteins that includes evaluating the quality of the crude and refined models. We used the ZLab server (https://zlab.umassmed.edu/bu/rama/) and the ProSA-web server (https://prosa.services.came.sbg.ac.at/prosa.php) to assess the quality of the initial and refined models. The Ramachandran plot is generated by the Zlab server and reveals the two-dimensional distribution of the torsion angles-phi (φ) and psi (ψ)-of the amino acids contained in the protein backbone [60,61]. The ProSA-web server computes an overall quality score for the protein structure. If this score falls outside of the normal range for native proteins, the structure most likely contains errors [62,63].

Identification of discontinuous B-cell epitopes
The ElliPro tool from the IEDB server (http://tools.iedb.org/ellipro/) was used to predict discontinuous B-cell epitopes in the vaccine's refined 3D model. ElliPro predicts and visualizes antibody epitopes in a given protein structure using Thornton's method, a residue clustering algorithm, the MODELLER program, and the Jmol viewer [64].

Disulfide engineering of the multi-epitope vaccine
Disulfide engineering is a method for forming disulfide bonds in protein structures. Disulfide bonds promote the stability of the folded protein structure by decreasing conformational entropy and increasing the free energy of the denatured state [65]. The Disulfide by Design 2.13 server (http://cptweb.cpt.wayne.edu/DbD2/) was used to identify residue pairs in the vaccine construct that had the potential to mutate to cysteine and form a disulfide bond [66].

Molecular docking of vaccine construct with TLR4
Molecular docking is an in silico method for predicting the binding affinity and orientation of a receptor and its ligand [67]. The Cluspro 2.0 server (https://cluspro.bu.edu/login.php) was used to perform molecular docking between the designed vaccine and Toll-Like receptor 4 (TLR4) (PDB ID: 4G8A). This server runs the following three computational steps to perform molecular docking: (1) rigid body docking is carried out using billions of conformations as samples, (2) clustering of the 1000 lowest energy structures using root-mean-square deviation (RMSD) to find the largest clusters that will represent the most likely models of the complex, and (3) chosen structures are refined using energy minimization [68][69][70][71]. The LigPlot software was used to generate schematic diagrams of bonds formed between vaccine construct residues and TLR4 [72].

Molecular dynamics simulation of the TLR4-vaccine complex
Molecular dynamics (MD) simulation is a powerful computational method for determining the motions of proteins at the atomic scale over time, based on a general model of the physics governing interatomic interactions [73]. In this study, we used the GROMACS 2019.6 software to run an MD simulation [74]. The selected complex from the molecular docking step was used as input for the MD simulation process. The input structure was prepared using the ff99SB force field. This structure was solved in a cubic box of TIP3P water molecules using the gmx solvate software, and then Na + and Cl − ions were introduced to the protein's surface to neutralize the total charge of the system. The system's energy was minimized using the steepest descent algorithm. The system's temperature was gradually increased from 0 to 300 K for 200 ps at constant volume, and the system was then brought to equilibrium at constant pressure. The MD simulation was run at 300 K and for 40 ns.

Calculation of binding free energy
The MMPBSA method and the gmx MMPBSA program were used to calculate the binding free energies for the vaccine, TLR4, and the vaccine-TLR4 complex. In this method, Poisson-Boltzmann (PB) and generalized Born (GB) models were used to calculate the binding free energy [75], and 1000 frames at regular intervals were extracted from the simulation trajectory and analyzed.

Codon optimization and in silico cloning of the multi-epitope vaccine
The codon optimization approach can improve the expression efficiency of foreign genes in the host [76]. The Java Codon Adaptation Tool (JCat) (http://www.jcat.de/) was used to perform back translation, codon optimization, and determine the vaccine sequence's codon adaptation index (CAI) value and GC content [77]. In this analysis, the protein sequence of the designed vaccine was submitted to this server as input, and E. coli (strain K12) was chosen as the host organism [46]. Three additional options were selected to avoid the rho-independent transcription termination, prokaryote ribosome binding site, and cleavage sites of restriction enzymes. Following that, the HindIII and BamHI restriction enzyme recognition sequences were introduced to the N-and C-terminal of the vaccine sequence, respectively. The sequence was then inserted into the expression pET28a (+) vector using the SnapGene 3.2. 1sofware.

Ethics statement
The ethical committee of Shahrekord University of Medical Sciences approved this study with the number: IR.SKUMS.REC.1401.065.

Retrieval of protein sequences and identification of conserved regions
A total of seven appropriate sequences for the envelope protein A28 homolog and eighteen suitable sequences for the cell surface-binding protein were selected (S1 and S2 Data). The Clustal Omega then performed MSA to identify the conserved regions of target proteins. The conserved regions were selected based on the protein sequence's lack of gaps and the highest number of similar amino acids.

Prediction and screening of T-cell and linear B-cell epitopes
Using the ProPred-I server, 78 CTL epitopes for the cell surface-binding protein and 56 CTL epitopes for the envelope protein A28 homolog were predicted. Among them, 12 CTL epitopes of the cell surface-binding protein and 21 CTL epitopes of the envelope protein A28 homolog that was capable of binding to at least three MHC class-I alleles were tested for antigenicity, toxicity, and allergenicity using the VaxiJen v2.0, ToxinPred, and AllerTOP v. 2.0 servers, respectively. The screened epitopes were then checked for their presence in the conserved regions of the target proteins. Finally, 2 CTL epitopes for the cell surface-binding protein (S1 Table) and 7 CTL epitopes for the envelope protein A28 homolog (S1 and S2 Tables) were chosen. Here, we predicted 38 HTL epitopes for the cell surface-binding protein and 25 HTL epitopes for the envelope protein A28 homolog using the ProPred server. Antigenicity, nonallergenicity, toxicity, presence in conserved regions, and IFN-γ and IL-4 production were assessed for 20 HTL epitopes of the cell surface-binding protein and 16 HTL epitopes of the envelope protein A28 homolog that could bind to at least three MHC class-II alleles. In total, 3 HTL epitopes were chosen for each target protein based on the mentioned parameters (S3 and S4 Tables). Moreover, the ABCpred server predicted 19 and 12 LBL epitopes for the cell surface-binding protein and envelope protein A28 homolog, respectively. After evaluation, 4 LBL epitopes for the cell surface-binding protein and 1 LBL epitope for the envelope protein A28 homolog with a score of 0.8 or higher were found to be antigenic, non-allergenic, and nontoxic in the conserved regions (S5 and S6 Tables).

Multi-epitope vaccine construction
Overlapping epitopes were removed to avoid epitope duplication in the vaccine construct's linear structure. The final vaccine structure consisted of 1 adjuvant (CTxB), 7 CTL, 4 HTL, and 5 LBL epitopes linked together with EAAAK, AAY, GGGGS, and KK linkers, in that order (Fig 2). The epitopes that constituted the vaccine's structure are listed in Table 1.

Population coverage analysis
To estimate population coverage of selected CTL and HTL epitopes, the IEDB population coverage analysis tool was employed. The CTL and HTL epitopes in the multi-epitope vaccine can cover 87.67% and 64.07% of the world population, respectively. The population coverage for combined CTL and HTL epitopes was estimated to be 95.57%. The population coverage of CTL and HTL epitopes when used in combination is highest in Europe (98.

Prediction of antigenicity, allergenicity, solubility, and physicochemical properties of the vaccine construct
The VaxiJen 2.0 server predicted the antigenicity of the designed vaccine to be 0.6355 with a virus model and 0.6107 with a bacteria model at a threshold of 0.4, while the ANTIGENpro server predicted the probability of antigenicity of the vaccine to be 0.774680. The vaccine construct was identified as non-allergen by the AllerTOP v.2.0 server. Furthermore, the Protein-Sol server predicted that our proposed vaccine was soluble, with a score of 0.455. The Prot-Param tool predicted the physiochemical characteristics of the vaccine candidate, which are displayed in Table 2. 3D structure modeling, refinement, and validation of the vaccine candidate I-TASSER server predicted the five 3D structure models of our vaccine construct using several threading templates (PDB Hit: 1LTR, 7P1A, 7LVB, 4L6T, 7QQK, 3KIG, 7Q0D, and 5ELB). The Z-score values ranged from 1 to 8.82, revealing that all of the threading templates were properly aligned; a Z-score greater than 1 indicates good alignment, and vice versa. Models 1-5 had C-scores of -1.29, -2.18, -2.45, -2.24, and -4.02, respectively. The C-score typically ranges from -5 to 2, with a higher C-score indicating a more confident model. As a result, we chose model 1 with a C-score of -1.29 for the refinement process (Fig 4). The GalaxyRefine server was used to refine the selected model, generating five refined models with different quality assessment parameters (Table 3). Higher GDT-HA and Rama favored values, as well as lower RMSD, MolProbity, Clash score, and Poor rotamers values, indicate that the models are of higher quality [78]. Model 2 was chosen for further evaluation based on the parameters  mentioned above (Fig 4). The Zlab server and the ProSA-web server were used to analyze the overall quality of the vaccine construct's three-dimensional structure before and after refining. The Ramachandran plot analysis of the initial model revealed that 83.188%, 11.304%, and 5.507% of residues were present in the highly preferred, preferred, and questionable regions,  respectively (Fig 5A), whereas, in the refined model, these values changed to 94.783%, 4.638%, and 0.580%, respectively (Fig 5B). The ProSA web server calculated the z-scores for the initial model to be -3.48 (Fig 5C), which changed to -3.67 after refinement (Fig 5D).

Identification of discontinuous B-cell epitopes
The ElliPro server predicted eight discontinuous B-cell epitopes with scores greater than 0.5 in the vaccine construct (Fig 6). The predicted epitopes ranged in size from 5 to 79 amino acids (Table 4).

Disulfide engineering of the multi-epitope vaccine
The Disulfide by Design 2.13 server identified 24 residue pairs in the vaccine construct's refined model that could form a disulfide bond ( Table 5). The χ3 peaks in 1505 native disulfide bonds in 331 non-homologous proteins have been observed at -87 and +97 degrees and approximately 90% of naturally formed disulfide bonds have an energy value of less than 2.2  kcal/mol [66]. Based on a χ3 angle between -87˚and +97˚and an energy value less than 2.2 kcal/mol, only three residue pairs were chosen for disulfide bond formation. These residue pairs were LYS 290-ALA 306, PHE 146-LYS 341, and SER 345-ALA 350 (Fig 7).

Molecular docking of vaccine construct with TLR4
Molecular docking between the vaccine and TLR4 was done using the Cluspro 2.0 server. This server generated 11 clusters and calculated the number of members and the corresponding lowest energy for each cluster. The comparison of all clusters revealed that cluster 1 has the strongest interaction energy (the most negative) with the lowest energy of -999.8 kcal/mol and 29 members (Fig 8). The LigPlot tool illustrated the bonds formed between the amino acids in the vaccine and TLR4 (Fig 9). The amino acids involved in hydrogen bonds, as well as their lengths, are listed in Tables 6 and 7. Comparing the quality of the initial and refined models of the vaccine construct's three-dimensional structure using the Zlab server and the ProSA-web server. (A) The Ramachandran plot analysis shows that in the initial model, 83.188%, 11.304%, and 5.507% of the residues are found in the highly preferred (green crosses), preferred (brown triangles), and questionable (red circles) regions, respectively; (B) whereas in the refined model, 94.783%, 4.638%, and 0.580% of the residues are found in the highly preferred, preferred, and questionable regions, respectively.

Molecular dynamics simulation of the TLR4-vaccine complex
The MD simulation results obtained from GROMACS 2019.6 software were analyzed as RMSD and RMSF using gmx rms and gmx rmsf modules, respectively. The RMSD parameter of the structures formed during MD simulation in the time dimension is an appropriate and  widely used measure for assessing the structural stability of the protein [79]. The RMSD plot of the TLR4 indicated values in the 0.16-0.52 nm range. The vaccine's RMSD plot showed an upward trend at the beginning of the simulation, peaking at 0.84 nm after about 33 ns and fluctuating slightly until the simulation ended ( Fig 10A). RMSF is a calculation of protein structure  residue fluctuation during an MD simulation. Because the two chains A and B from the TLR4 have the same sequence, the RMSF plot of these two chains is shown to more accurately study the effect of the vaccine construct binding on the flexibility of chain A. The residues at both ends of chain A have shown less flexibility than chain B, while the flexibility of other regions in the two chains is the same and does not show a noticeable change. The vaccine's RMSF plot revealed that the majority of the vaccine's residues were highly flexible (Fig 10B).

Calculation of binding free energy
The binding free energies of TLR4, vaccine, and TLR4-vaccine complex were calculated using the MM-GBSA and MM-PBSA methods. In the MM-PBSA analysis, the total binding free energy was -53718.65 kcal/mol for TLR4, -33558.07 kcal/mol for the vaccine, and -87411.4 kcal/mol for the TLR4-vaccine complex. Moreover, the total binding free energy in the MM-GBSA was -54463.47 kcal/mol for TLR4, -34327.63 kcal/mol for the vaccine, and -88867.65 kcal/mol for the TLR4-vaccine complex. According to the calculated values, the gas phase energy contribution was equal and substantial in both methods. The contribution of each energy component is shown in Table 8.

Codon optimization and in silico cloning of the multi-epitope vaccine
The vaccine protein sequence was reverse translated into a nucleotide sequence with a length of 1128 bp using the JCat server. This server estimated the optimized sequence's GC content to be 45.30% and its CIA value to be 0.97. The vaccine nucleotide sequence was then in silico cloned into the multiple cloning site (MCS) of pET28a (+) between HindIII (173) and BamHI (1307) restriction sites, resulting in a recombinant plasmid with a length of 6478 bp (Fig 11).

Discussion
Despite the fact that the MPXV is endemic in West and Central African countries, outbreaks in other countries with no known epidemiological links to West or Central Africa have occurred [80]. Because the MPXV is so closely related to smallpox, the smallpox vaccine is at least 85% effective in preventing MPX, but a vaccine that specifically targets the MPXV is needed [81]. The multi-epitope vaccines have a high potential to elicit humoral and cellular immune responses; these vaccines are composed of epitopes rather than large proteins or whole genomes, protecting the host from excessive antigenic load and allergic reactions [82].
In recent years, a large number of multi-epitope vaccines against various pathogens have been designed, including human papillomavirus (HPVs) [83], Zika virus [84,85], Ebola virus [86,87], SARS-CoV-2 [88][89][90][91][92], Klebsiella pneumonia [93], Helicobacter pylori [79,94], Staphylococcus aureus [95], Leishmania [96,97], and Cystic echinococcosis [98]. Several studies on multi-epitope vaccine design against MPXV have been carried out so far. In the study conducted by Aiman et al. (2022), selected epitopes from MPXVgp167, MPXVgp028, and MPXVgp105 proteins along with 4 different adjuvants (HBHA protein, βdefensins, 50S ribosomal protein L7/L12 adjuvants, and HBHA conserved peptide sequences) were included in four vaccine constructs [99]. In the study of Ullah et al. (2022), to design a multi-epitope vaccine against MPXV, three extracellular proteins, including a cupin domaincontaining protein, ABC transporter ATP-binding protein, and DUF192 domain-containing protein were selected as epitope prediction targets, as well as the CTxB was used as an adjuvant in the vaccine construct [100]. The multi-epitope vaccine designed by Shantier et al. (2022) included T-cell and B-cell epitopes selected from the cell surface-binding proteins, as well as Table 6. List of residues involved in the formation of hydrogen bonds between TLR4 (chains A) and the vaccine.  Table 7. List of residues involved in the formation of hydrogen bonds between TLR4 (chain C) and the vaccine.   protein L7/L12, b-defensin, LT-IIC, CTB, RS09) resulted in five different vaccine constructs [104]. In the study of Jin, et al. (2023) a vaccine construct was designed that included B-and T-cell epitopes from B13R, Abundant component of virosome, Bifunctional zinc finger-like protein/E, Ankyrin, DNA-directed RNA polymerase subunit, A18L, Bifunctional hemagglutinin/type-I membrane, IMV membrane protein, Interleukin-1 receptor antagonist-like protein, and CC-type chemokine binding protein. They also appended a beta-defensin-2 sequence as an adjuvant to the vaccine construct's N-terminal [105]. In comparison to published articles, our study is unique in that it combined the epitopes of the cell surface-binding protein and envelope protein A28 homolog for the first time with CTxB (adjuvant) to construct a multiepitope vaccine against MPXV. Because the cell surface-binding protein and envelope protein A28 homolog are important for virus infectivity, we considered them for epitope prediction [106]. Epitope mapping of the aforementioned proteins was performed, and epitopes interacting with at least three MHC class-I and MHC class-II alleles were chosen for further screening. We planned to design a universal multi-epitope vaccine against the MPXV. For this purpose, in addition to evaluating epitopes for parameters such as antigenicity, toxicity, and allergenicity, as well as induction of IFN-γ and IL-4 (only for HTL epitopes), we also considered their presence in conserved regions from all virus strains available in NCBI. The multi-epitope vaccine was constructed by assembling the selected BCL, HTL, and CTL epitopes and appropriate linkers as well as the CTxB (as adjuvant). The CTxB is a non-toxic component of cholera toxin that has a high affinity for the monosialotetrahexosylganglioside (GM1) receptor, which is found in many cell types, including gut epithelial cells, antigen-presenting cells, macrophages, dendritic cells, and B cells [107]. As a result, CTxB has been widely employed as an adjuvant in vaccine design to boost immune responses [108][109][110]. In this study, like many other studies, the AAY linker was used to connect CTL epitopes [85,111,112], the GGGGS linker to connect HTL epitopes [113,114], the KK linker to connect LBL epitopes [115][116][117], and the EAAAK linker to connect adjuvant [116,118,119]. The AAY linkers help to increase epitope presentation while decreasing the formation of junctional epitopes [120]. The GGGGS linker, which is made up of small or polar amino acids like Gly and Ser, provides good flexibility and solubility and is an excellent choice for fusion protein domains that require specific movements or interactions [121]. The KK linker is a target for cathepsin B, a key protease in antigen processing. This linker helps in the reduction of junctional immunogenicity by preventing antibody induction for the peptide sequence formed when the individual epitopes are linearly connected [120]. The EAAAK linker is a rigid alpha-helix-forming peptide linker that functions as a domain spacer in fusion proteins [122].

TLR4 (C) Vaccine (E) Bond Length (Å)
According to population coverage analysis, the vaccine construct covers 95.57% of the worldwide population, while the vaccine designed in study of Waqas et al. covers 93.62% of the worldwide population [102]. The high antigenicity scores predicted by the VaxiJen 2.0 server and the ANTIGENpro server showed that the proposed vaccine is effective enough to elicit strong immune responses in the body. Furthermore, the vaccine construct's non-allergic behavior was validated using the AllerTOP v. 2.0 server. The vaccine has good solubility, according to the Protein-Sol server's solubility score (0.455), which was higher than the threshold (0.45). The solubility score of the vaccine designed in the study of Shantier et al. (2022) was lower than the threshold, so in terms of solubility, the vaccine designed in the present study is preferable to this vaccine [101]. Evaluating the physicochemical properties of the vaccine structure provides us with relative knowledge of the vaccine's physicochemical properties, which facilitates in vitro and in vivo vaccine evaluations. The vaccine construct had a molecular weight of 41.97 kDa, which is ideal because proteins with a molecular weight of less than 110 kDa are easier and faster to express and purify than heavier proteins [95]. The proposed vaccine's theoretical pI was determined to be 9.48, implying its alkaline nature. Our vaccine candidate's half-life was estimated to be 30 hours in mammalian reticulocytes, more than 20 hours in yeast, and more than 10 hours in E.coli, indicating that the vaccine is exposed to the immune system for an extended period of time [123]. In terms of half-life, our vaccine outperformed the vaccine designed by Waqas et al., which had a half-life of 7.2 hours in mammalian reticulocytes (in vitro) [102]. The half-life of our vaccine matched the half-life of the vaccine construct designed in the study by Ullah et al. [100] and Shantier et al. [101]. The vaccine instability index was determined to be 38.70, which is less than 40, indicating that the vaccine candidate is stable under standard conditions [123]. The vaccine's aliphatic index was computed to be 91.12, indicating that the vaccine is stable over a wide temperature range [124]. The aliphatic index value of this vaccine was higher than the aliphatic index of the vaccine designed in studies by Waqas et al. (74.63) [102], Zaib et al. (79.90) [103], and Shantier et al. (65.75) [101]. The GRAVY score of the multi-epitope vaccine was calculated to be -0.044; this parameter indicates hydrophilicity and is related to protein solubility, and its negative value indicates that the vaccine interacts better with water molecules [32,125].
After modeling, the vaccine's 3D structure was refined to optimize its quality and bring it closer to experimental accuracy. The Ramachandran plot revealed that in the initial model, 83.188% of the residues were present in the highly preferred region, while after refinement, the number of residues in this region increased to 94.783%. The z-scores for the initial and refined models were calculated by the ProSA web server as -3.48 and -3.76, respectively. A higher negative value of this parameter indicates that the structure is of higher quality [22]. Humoral immunity along with cellular immunity, is critical in fighting pathogens in the body [126]. The B-cells which produce antibodies are responsible for the body's humoral immunity [127]. Discontinuous B-cell epitopes play a key role in the antigen-antibody response [128]. Our vaccine candidate has the potential to induce robust humoral immunity because it has a large number of discontinuous B-cell epitopes.
Various Toll-like receptors (TLRs) are involved in the initial interaction of host cells with invading viruses, regulate virus replication and host responses, and ultimately affect virus pathogenesis [129]. Hutchens et al. revealed that in vaccinia-infected mice, TLR4 recognizes a viral ligand rather than an endogenous ligand, and mice with a nonfunctional mutant version of TLR4 had higher viral replication, hypothermia, and mortality than control animals. The findings of this work demonstrated that TLR4 promotes a protective innate immune response against the vaccinia virus, which can help to develop new vaccines and therapeutic agents against poxviruses [130], therefore, we used TLR4 as a receptor in molecular docking analysis. The results of the molecular docking show that our vaccine has a strong affinity for TLR4 and can stimulate both innate and adaptive immunity. The RMSD plot revealed the vaccine-TLR4 complex's stability during the MD simulation The vaccine's binding to chain A of TLR4 has reduced chain A's flexibility in comparison to chain B, as shown by the RMSF plot. Furthermore, the high flexibility of the vaccine's amino acids 40-100 in comparison to its other amino acids can be attributed to the fact that this region has no interaction with the receptor protein and can move freely. The computation of the binding free energies confirmed that the TLR4-vaccine complex is more stable than each of its parts individually.
The CAI value and GC content of our vaccine were calculated by the JCat server to be 0.97 and 45.30%, respectively. CAI value between 0.8 and 1 [131], as well as GC content between 30% and 70% [128], are considered optimal for gene expression in the target organism. Although the findings of this study are encouraging, it is still essential to test the vaccine candidate both in vitro and in vivo.

Conclusion
Following recent WHO reports of confirmed MPXV cases in non-endemic regions, global concern about the possibility of a new pandemic has been raised. There is still no licensed vaccine specifically for MPXV. Because computational methods allow researchers to develop an effective and safe vaccine in less time and at a lower cost, we used this strategy to develop a multi-epitope vaccine against MPXV in this study. We believe our vaccine candidate can be a universal vaccine capable of inducing cellular and humoral immunity because we incorporated the best T-cell and B-cell epitopes from the conserved regions of the target proteins into the vaccine construct. Furthermore, the proposed vaccine's high affinity for the innate immune receptor (TLR4) indicates that it could be able to stimulate both innate and adaptive immunity against pathogen infection. The findings of this study are encouraging, but they are based on computational methods that will never be able to replace laboratory validation; thus, we recommend that this designed vaccine construct be tested in animal models in the future.